home *** CD-ROM | disk | FTP | other *** search
/ TPUG - Toronto PET Users Group / TPUG Users Group CD / TPUG Users Group CD.iso / PET / S-Super PET / (s)t2.d64 / PERIODOGRAM.FTN < prev    next >
Text File  |  2009-01-18  |  8KB  |  304 lines
  1. *           Program Periodogram.ftn
  2. *
  3. *  This program was typed from the book:
  4. *
  5. *    FOURIER ANALYSIS OF TIME SERIES: AN INTRODUCTION
  6. *    by Peter Bloomfield
  7. *    John Wiley & Sons, 1976
  8. *
  9. *   This program computes the discrete Fourier transform
  10. *   of a series,  and the periodogram.
  11. *
  12. *
  13.       real x(1024) , y(1024)
  14.       character fnamein,realpartfile,imagpartfile,periodogramfile
  15.       maxsize=1024
  16.  
  17.       do 10 i = 1 , maxsize
  18.          x(i) = 0.0
  19.          y(i) = 0.0
  20.    10 continue
  21.  
  22.       write(6,1)
  23.     1 format("1 Enter input filename (up to 16 characters)")
  24.       read, fnamein
  25.       print, "Output file name for the periodogram"
  26.       read, periodogramfile
  27.  
  28.       call datin (x , n , start , step , 7 , fnamein)
  29.  
  30.       if n .gt. maxsize
  31.           print,"Number of data points is ",n," which exceeds maximum size"
  32.           stop
  33.       endif
  34.  
  35.       call detrnd( x, n, 0 )
  36.  
  37.       np2 = 1
  38.    20 np2 = np2*2
  39.       if (np2 .lt. n) go to 20
  40.  
  41.       inv = 0
  42.       call ft01a (x , y , np2 , inv)
  43.  
  44.       if (inv .eq. 0) go to 30
  45.          write(6,7)
  46.     7    format("0    TROUBLE IN ft01a  --  np2 OUT OF RANGE")
  47.          stop
  48.    30 continue
  49.  
  50.       npgm = np2/2 + 1
  51.       con = (2.0*float(np2)/float(n))**2
  52.       do 40 i = 1 , npgm
  53.    40 x(i) = (x(i)**2 + y(i)**2)*con
  54.  
  55.       call datout (x , npgm , 0.0 , 1.0 , 9 , periodogramfile)
  56.       stop
  57.       end
  58.  
  59.  
  60.       subroutine datin ( x , n , start , step , m , fname)
  61. *
  62. *  This subroutine is used to input a series of values
  63. *  (in run time format) and some associated quantities
  64. *  (in fixed format).  The first four data records are -
  65. *
  66. *  1   A header record (72 bytes)
  67. *  2   Value of n (i5)
  68. *  3   The data format (72 bytes)
  69. *  4   Start and step  (2f10.5)
  70. *
  71. *  Parameters are -
  72. *
  73. *  Name      Type                     Value
  74. *                         On Entry                On Return
  75. *
  76. *  x      Real array      Not used              The series
  77. *
  78. *  n      Integer         Series length         Unchanged
  79. *
  80. *  start  Real            Not used              Time value at the
  81. *                                               first data point
  82. *
  83. *  step   Real            Not used              Time increment
  84. *                                               between data points
  85. *
  86. *  m      Integer         Logical unit number   Unchanged
  87. *
  88. *  fname  Character       input file name      Unchanged
  89. *
  90.       real x(1024)
  91.       character head , form , fname
  92.       open(unit=m,file=fname,access="sequential",recl=80)
  93.       read(m,1) head,n,form,start,step
  94.     1 format(a72/i5/a72/2f10.5)
  95.       write(6,2) head,n,form,start,step
  96.     2 format("0The data header reads -"/1x,a72/
  97. &            " The series length is",i6/
  98. &            " The data format is -"/1x,a72/
  99. &            " time origin is",f11.5,
  100. &            ",  time increment is",f11.5)
  101.       read(m,form) (x(i),i=1,n)
  102.       close(unit=m)
  103.       return
  104.       end
  105.  
  106.       subroutine detrnd (x , n , k)
  107. *
  108. *  This subroutine computes and subtracts from the
  109. *  series  x  either the series mean or the least squares
  110. *  straight line.  Note that these are the least squares
  111. *  polynomials of degree  0  and  1,  respectively.
  112. *  Parameters are -
  113. *
  114. *  Name     Type                      Value
  115. *                         On Entry              On Return
  116. *
  117. *  x     Real array    The series             Detrended series
  118. *
  119. *  n     Integer       Series length          Unchanged
  120. *
  121. *  k     Integer       Degree of polynomial   Unchanged
  122. *                      to be fitted
  123. *
  124.       real x(n)
  125.       sumx = 0.0
  126.       do 20 i = 1 , n
  127.          sumx = sumx + x(i)
  128.    20 continue
  129.       xbar = sumx/float(n)
  130.       do 30 i = 1 , n
  131.          x(i) = x(i) -xbar
  132.    30 continue
  133.       if (k .le. 0) return
  134.       tbar  = float(n+1)/2.0
  135.       fn    = n
  136.       sumtt = fn*(fn*fn-1.0)/12.0
  137.       sumtx = 0
  138.       do 40 i = 1 , n
  139.          sumtx = sumtx + x(i)*(float(i)-tbar)
  140.    40 continue
  141.       beta = sumtx/sumtt
  142.       do 50 i = 1 , n
  143.          x(i) = x(i) - beta*(float(i)-tbar)
  144.    50 continue
  145.       return
  146.       end
  147.  
  148.  
  149.       subroutine datout (x , n , start , step , m , fname)
  150. *
  151. *  This subroutine outputs the series  x  in the format
  152. *  expected by subroutine  datin.   The header record is
  153. *  blank except for sequencing.   Parameters are -
  154. *
  155. *  Name       Type           Value (none are changed)
  156. *
  157. *  x        Real Array       The series
  158. *
  159. *  n        Integer          Series length
  160. *
  161. *  start    Real             Time origin
  162. *
  163. *  step     Real             Time increment
  164. *
  165. *  m        Integer          Logical unit number
  166. *
  167. *  fname    Character        Output file name
  168. *
  169.       real x(n),z(4)
  170.       character fname
  171.       open(unit=m,file=fname,access="sequential",recl=80)
  172.       do 7 i = 1 , 4
  173.     7 z(i) = 0.0
  174.       k = 1
  175.       write(m,1) fname , k
  176.     1 format(1x,a74,i5)
  177.       k = 2
  178.       write(m,2) n,k
  179.     2 format(i5,70x,i5)
  180.       k = 3
  181.       write(m,3) k
  182.     3 format(" (5e15.8)",66x,i5)
  183.       k = 4
  184.       write(m,4) start,step,k
  185.     4 format(2f10.5,55x,i5)
  186.       l = 5
  187.       ihi = 0
  188.    10 ilo = ihi + 1
  189.       ihi = ihi + l
  190.       k   = k + 1
  191.       if (ihi .gt. n) go to 20
  192.       write(m,5) (x(i),i=ilo,ihi),k
  193.     5 format(5e15.8,i5)
  194.       go to 10
  195.    20 if (ilo .gt. n) go to 30
  196.       lim = ihi - n
  197.       write(m,5) (x(i),i=ilo,n),(z(i),i=1,lim),k
  198.    30 close(unit=m)
  199.       return
  200.       end
  201.  
  202.       subroutine ft01a (xr,xi,n,inv)
  203. *
  204. *    This subroutine implements the Sande-Tukey radix-2
  205. *    Fast Fourier Transform.  Either the direct or
  206. *    the inverse transform may be computed.  Parameters are
  207. *
  208. *  Name     Type                         Value
  209. *                          On Entry                 On Return
  210. *
  211. *  xr    real array   real part of the            real part of the
  212. *                     series                      transform
  213. *
  214. *  xi    real array   imaginary part of           imaginary part
  215. *                     the series                  of the transform
  216. *
  217. *  n     integer      the series length           unchanged
  218. *
  219. *  inv   integer      0  for direct transform     inv is set to
  220. *                     1  for inverse transform    -1 to indicate
  221. *                                                 error return
  222. *
  223. *   The direct transform is -
  224. *
  225. *             n
  226. *      (1/n) sum x(t)*exp(-2*pi*i*(t-1)*(j-1)/n),   j = 1 to n
  227. *            t=1
  228. *
  229. *   The inverse transform is -
  230. *
  231. *             n
  232. *            sum x(t)*exp( 2*pi*i*(t-1)*(j-1)/n),   j = 1 to n
  233. *            t=1
  234. *
  235.       real xr(n),xi(n),ur(15),ui(15)
  236.  
  237.       ur(1) = 0.0
  238.       ui(1) = 1.0
  239.       do 110 i = 2 , 15
  240.          ur(i) = sqrt((1.0+ur(i-1))/2.0)
  241.  110  ui(i) = ui(i-1)/(2.0*ur(i))
  242.  
  243.  120  if (n .gt. 0  .and.  n .le. 2**14) go to 130
  244.       inv = -1
  245.       return
  246.  
  247.  130  n0 = 1
  248.       ii = 0
  249.  140  n0 = n0 + n0
  250.       ii = ii + 1
  251.       if (n0 .lt. n) go to 140
  252.       i1 = n0/2
  253.       i3 = 1
  254.       i0 = ii
  255.  
  256.       do 260 i4 = 1 , ii
  257.          do 250 k = 1, i1
  258.             wr = 1.0
  259.             wi = 0.0
  260.             kk = k - 1
  261.             do 230 i = 1 , i0
  262.                if (kk .eq. 0) go to 240
  263.                if (mod(kk,2) .eq. 0) go to 230
  264.                j0 = i0 - i
  265.                ws = wr*ur(j0) - wi*ui(j0)
  266.                wi = wr*ui(j0) + wi*ur(j0)
  267.                wr = ws
  268.  230        kk = kk/2
  269.  240        if ( inv .eq. 0 ) wi = -wi
  270.             l = k
  271.             do 250 j = 1 , i3
  272.                l1 = l + i1
  273.                zr = xr(l) + xr(l1)
  274.                zi = xi(l) + xi(l1)
  275.                z      = wr*(xr(l)-xr(l1))-wi*(xi(l)-xi(l1))
  276.                xi(l1) = wr*(xi(l)-xi(l1))+wi*(xr(l)-xr(l1))
  277.                xr(l1) = z
  278.                xr(l)  = zr
  279.                xi(l)  = zi
  280.  250     l = l1 + i1
  281.          i0 = i0 - 1
  282.          i3 = i3 + i3
  283.  260  i1 = i1/2
  284.  
  285.       um = 1.0
  286.       if ( inv .eq. 0 ) um = 1.0/float(n0)
  287.       do 310 j = 1 , n0
  288.          k  = 0
  289.          j1 = j - 1
  290.          do 320 i = 1 , ii
  291.             k = 2*k + mod(j1,2)
  292.  320     j1 = j1/2
  293.          k = k + 1
  294.          if ( k .lt. j ) go to 310
  295.          zr    = xr(j)
  296.          zi    = xi(j)
  297.          xr(j) = xr(k)*um
  298.          xi(j) = xi(k)*um
  299.          xr(k) = zr*um
  300.          xi(k) = zi*um
  301.  310  continue
  302.       return
  303.       end
  304.